---
title: 'Section 4: The 2019 Istanbul Mayoral Election Quantitative Analysis'
author: "Milan Svolik"
date: "4/20/2023"
output:
  pdf_document: default
  html_document:
    df_print: paged
---

```{r setup, include=FALSE}
rm(list = ls(all = TRUE))
library(tidyverse)
library(stargazer)
library(lmtest)
library(reshape2)
library(foreign)
library(ggplot2)
library(broom)
library(estimatr)

set.seed(5877)
```

## Load the data
```{R}
load("data/df_2019istanbul_results.RData")
```

### Prepare the quantities that will be plotted
```{R}
# AGGREGATE THE RESULTS AT THE NEIGHBORHOOD LEVEL
df_results_neighbourhood <- df_results %>% 
  group_by(district_id, district_name, neighborhood_id, neighborhood_name) %>% 
  summarize(AKP_march=sum(AKP_march), CHP_march=sum(CHP_march), ABS_march=sum(ABS_march), AKP_june=sum(AKP_june), CHP_june=sum(CHP_june), ABS_june=sum(ABS_june))

# AKP VOTE, CHP VOTE, ABSTENTIONS AS A SHARE OF REGISTERED VOTERS
df_results_neighbourhood <- df_results_neighbourhood %>%
  mutate(
    REG_march = AKP_march + CHP_march + ABS_march,  
    AKP_share_march = AKP_march/REG_march, 
    CHP_share_march = CHP_march/REG_march,
    ABS_share_march = ABS_march/REG_march, 
    REG_june = AKP_june + CHP_june + ABS_june,  
    AKP_share_june = AKP_june/REG_june, 
    CHP_share_june = CHP_june/REG_june,
    ABS_share_june = ABS_june/REG_june, 
  )

# SHIFTS IN AKP VOTE SHARE, CHP VOTE SHARE, AND ABSTENTION RATE
df_results_neighbourhood <- df_results_neighbourhood %>%
  mutate(
    AKP_share_diff = AKP_share_june - AKP_share_march,
    CHP_share_diff = CHP_share_june - CHP_share_march, 
    ABS_share_diff = ABS_share_june - ABS_share_march
  )
```
# Figure 4.4: 
* Neighborhood-level distribution of March 2019 the AKP vote, CHP vote, and abstention as a fraction of registered voters
```{R}
df_results_neighbourhood %>%
  ggplot(aes(x = AKP_share_march)) +
  geom_histogram(mapping = aes(y=..count../sum(..count..)), binwidth = .05, colour = "black", fill = "white") +
  theme_minimal() + 
  xlim(0,1) + 
  labs(
       x = "March AKP vote as a fraction of registered voters", y = "Fraction",
       title = ""
    )
ggsave("figures, appendix/akp_march_histogram.pdf", width = 10, height = 6, units = "in", dpi = 600)

df_results_neighbourhood %>%
  ggplot(aes(x = CHP_share_march)) +
  geom_histogram(mapping = aes(y=..count../sum(..count..)), binwidth = .05, colour = "black", fill = "white") +
  theme_minimal() + 
  xlim(0,1) + 
  labs(
       x = "March CHP vote as a fraction of registered voters", y = "Fraction",
       title = ""
    )
ggsave("figures, appendix/chp_march_histogram.pdf", width = 10, height = 6, units = "in", dpi = 600)

df_results_neighbourhood %>%
  ggplot(aes(x = ABS_share_march)) +
  geom_histogram(mapping = aes(y=..count../sum(..count..)), binwidth = .05, colour = "black", fill = "white") +
  theme_minimal() + 
  xlim(0,1) + 
  labs(
       x = "March abstentions as a fraction of registered voters", y = "Fraction",
       title = ""
    )
ggsave("figures, appendix/abs_march_histogram.pdf", width = 10, height = 6, units = "in", dpi = 600)
```

# Figure 4.5
*Neighborhood-level distribution of shifts between March and June 2019 in the AKP vote, CHP vote, and abstention as a fraction of registered voters
```{R}
df_results_neighbourhood %>%
  ggplot(aes(x = AKP_share_diff)) +
  geom_histogram(mapping = aes(y=..count../sum(..count..)), binwidth = .01, colour = "black", fill = "white") +
  theme_minimal() + 
  xlim(-.2,.2) + 
  labs(
       x = "June-March shift in AKP vote", y = "Fraction",
       title = ""
    )
ggsave("figures, appendix/akp_diff_histogram.pdf", width = 10, height = 6, units = "in", dpi = 600)

df_results_neighbourhood %>%
  ggplot(aes(x = CHP_share_diff)) +
  geom_histogram(mapping = aes(y=..count../sum(..count..)), binwidth = .01, colour = "black", fill = "white") +
  theme_minimal() + 
  xlim(-.2,.2) + 
  labs(
       x = "June-March shift in CHP vote", y = "Fraction",
       title = ""
    )
ggsave("figures, appendix/chp_diff_histogram.pdf", width = 10, height = 6, units = "in", dpi = 600)

df_results_neighbourhood %>%
  ggplot(aes(x = ABS_share_diff)) +
  geom_histogram(mapping = aes(y=..count../sum(..count..)), binwidth = .01, colour = "black", fill = "white") +
  theme_minimal() + 
  xlim(-.2,.2) + 
  labs(
       x = "June-March shift in abstentions", y = "Fraction",
       title = ""
    )
ggsave("figures, appendix/abs_diff_histogram.pdf", width = 10, height = 6, units = "in", dpi = 600)
```
# Figure 4.6
*Neighborhood-level distribution of shifts between March and June 2019 in the AKP vote, CHP vote, and abstention as a fraction of registered voters
```{R}
df_results_neighbourhood %>%
  ggplot(aes(AKP_share_march, AKP_share_june)) +
  geom_abline(slope=1, col="gray50") +
  geom_point(col="blue")  + 
  geom_smooth(method=lm, se = FALSE, col="black") +
  theme_minimal() + 
  theme(aspect.ratio = 1) +
  coord_fixed(ratio=1, xlim = c(0,1), ylim = c(0,1)) +
  scale_x_continuous(name = "March 2019") +
  scale_y_continuous(name = "June 2019") +
  labs(title = "AKP Vote/Registered")

ggsave("figures, appendix/AKP, 45 degree.pdf", width = 6, height = 6, units = "in")

df_results_neighbourhood %>%
  ggplot(aes(CHP_share_march, CHP_share_june)) +
  geom_abline(slope=1, col="gray50") +
  geom_point(col="red")  + 
  geom_smooth(method=lm, se = FALSE, col="black") +
  theme_minimal() + 
  theme(aspect.ratio = 1) +
  coord_fixed(ratio=1, xlim = c(0,1), ylim = c(0,1)) +
  scale_x_continuous(name = "March 2019") +
  scale_y_continuous(name = "June 2019") +
  labs(title = "CHP Vote/Registered")

ggsave("figures, appendix/CHP, 45 degree.pdf", width = 6, height = 6, units = "in")

df_results_neighbourhood %>%
  ggplot(aes(ABS_share_march, ABS_share_june)) +
  geom_abline(slope=1, col="gray50") +
  geom_point(col="darkgreen")  + 
  geom_smooth(method=lm, se = FALSE, col="black") +
  theme_minimal() + 
  theme(aspect.ratio = 1) +
  coord_fixed(ratio=1, xlim = c(0,1), ylim = c(0,1)) +
  scale_x_continuous(name = "March 2019") +
  scale_y_continuous(name = "June 2019") +
  labs(title = "Abstention/Registered")

ggsave("figures, appendix/ABS, 45 degree.pdf", width = 6, height = 6, units = "in")
```
# Figure 4.7
*The distribution of neighborhood-level outcome shifts that point to vote-switching (top), backlash (middle), and disengagement (bottom)

### Prepare the data
```{R}
rm(list = ls(all = TRUE))
load("data/df_2019istanbul_results.RData")

## Aggregate results at the neighborhood level
df_results <- df_results %>% 
  mutate(
    ABS_march = ABS_march + THIRD_march,
    ABS_june = ABS_june + THIRD_june
    ) %>%
  select(-THIRD_march, -THIRD_june)

# AGGREGATE THE RESULTS AT THE NEIGHBORHOOD LEVEL
df_results_neighbourhood <- df_results %>% 
  group_by(district_id, district_name, neighborhood_id, neighborhood_name) %>% 
  summarize(
    AKP_march=sum(AKP_march), CHP_march=sum(CHP_march), ABS_march=sum(ABS_march), 
    AKP_june=sum(AKP_june), CHP_june=sum(CHP_june), ABS_june=sum(ABS_june), 
    REG_march = AKP_march + CHP_march + ABS_march,  
    AKP_share_march = AKP_march/REG_march, 
    CHP_share_march = CHP_march/REG_march,
    ABS_share_march = ABS_march/REG_march, 
    REG_june = AKP_june + CHP_june + ABS_june,  
    AKP_share_june = AKP_june/REG_june, 
    CHP_share_june = CHP_june/REG_june,
    ABS_share_june = ABS_june/REG_june, 
    AKP_share_diff = AKP_share_june - AKP_share_march,
    CHP_share_diff = CHP_share_june - CHP_share_march, 
    ABS_share_diff = ABS_share_june - ABS_share_march
  )
```

### SIMULATE THE NULL HYPOTHESIS
* Take the March fraction of citizens that voted for the AKP, the CHP, and abstained as the respective outcome probabilities under the null hypothesis
* Draw a large number of draws from the multinomial distribution for each neighborhood
* This simulation accounts for: 
  + differences in neighborhood size: the same percentage shift in an outcome is less likely to occur by chance in neighborhoods with a larger number of registered voters
  + the fact that a negative correlation between shifts in any two outcomes arises mechanically because the number of registered voters is constant across the two elections
```{r}
# A FUNCTION THAT RETURNS THE PERCENTILE OF THE VALUE x IN THE EMPIRICAL DISTRIBUTION OF sample
centile <- function(sample, x) {
              ecdf(sample)(x)  
}

# DATA LENGTH
neighborhoods_length <- nrow(df_results_neighbourhood)

# NUMBER OF DRAWS PER NEIGHBORHOOD
rep <- 10000

# LISTS TO SAVE SIMS INTO
vote_centiles_list <- vector("list", neighborhoods_length)

# LOOP OVER ALL NEIGHBORHOODS
for (i in 1:neighborhoods_length) {
  
  # ISOLATE MARCH AND JUNE OUTCOMES FOR NEIGHBORHOOD i
  df_neighborhood_i <- df_results_neighbourhood[i,]
  march_i <- c(df_neighborhood_i$AKP_march, df_neighborhood_i$CHP_march, df_neighborhood_i$ABS_march)
  june_i <- c(df_neighborhood_i$AKP_june, df_neighborhood_i$CHP_june, df_neighborhood_i$ABS_june)
  
  # THE NULL/MARCH PROBABILITIES 
  P_march_i <- march_i/sum(march_i)
  # rep DRAWS FROM THE NULL DISTRIBUTION
  NULL_draw_i <- rmultinom(n=rep, size=sum(june_i), prob=P_march_i)    

  # THE QUANTILE OF THE OBSERVED VALUE: HOW MUCH OF AN OUTLIER IS THE JUNE OUTCOME BY THE STANDARD OF THE MARCH BASED NULL_draw_i
  AKP_centile <- centile(NULL_draw_i[1,], june_i[1])
  CHP_centile <- centile(NULL_draw_i[2,], june_i[2])
  ABS_centile <- centile(NULL_draw_i[3,], june_i[3])

  # SAVE THE CENTILES FOR NEIGHBORHOOD i
  vote_centiles_list[[i]] <- c(AKP_centile, CHP_centile, ABS_centile)
  
  }
```

### Which June election outcomes depart from March ones too much to have occurred solely due to the idiosyncratic randomness inherent in elections? 
* identify all neighborhoods where at least one of the three June outcomes falls below the 2.5th or above the 97.5th percentile of the simulated draws
```{R}
# MATRIX WITH THE CENTILE COLUMNS: AKP, CHP, ABS
centiles <- matrix(unlist(vote_centiles_list), ncol=3, byrow=TRUE)

# CLASSIFY INTO "above", "below" (the 2.5th or above the 97.5th percentile of the simulated draws) or "neither"
outside <- ifelse(centiles<.025, "below", centiles)
outside <- ifelse(centiles>.975, "above", outside)
outside <- ifelse(!(centiles<.025 | centiles>.975) , "neither", outside)
```

### Construct the plots
```{r}
# CREATE A DUMMY FOR NEIGHBORHOODS THAT POINT TO VOTE-SWITCHING
akp_chp_sign <- ifelse(outside[,1]=="below" & outside[,2]=="above", 1, 0)
table(akp_chp_sign)
akp_chp_sign <- as.factor(akp_chp_sign)

# CREATE A DUMMY FOR NEIGHBORHOODS THAT POINT TO DISENCHANTMENT
akp_abs_sign <- ifelse(outside[,1]=="below" & outside[,3]=="above", 1, 0)
table(akp_abs_sign)
prop.table(table(akp_abs_sign))
akp_abs_sign <- as.factor(akp_abs_sign)
  
# CREATE A DUMMY FOR NEIGHBORHOODS THAT POINT TO BACKLASH
chp_abs_sign <- ifelse(outside[,2]=="above" & outside[,3]=="below", 1, 0)
table(chp_abs_sign)
prop.table(table(chp_abs_sign))
chp_abs_sign <- as.factor(chp_abs_sign)

df_plot <- df_results_neighbourhood %>%
  select(AKP_share_diff, CHP_share_diff, ABS_share_diff) %>%
  ungroup() %>%
  mutate(akp_chp_sign, akp_abs_sign, chp_abs_sign)

x_lim <- .25

df_plot %>%
  ggplot(aes(x=AKP_share_diff, y=CHP_share_diff, group=akp_chp_sign)) +
  geom_point(aes(shape=akp_chp_sign, col=akp_chp_sign, shape=akp_chp_sign)) +
  geom_hline(yintercept=0, size=1) +
  geom_vline(xintercept=0, size=1) +
  theme_minimal() + 
  coord_cartesian(xlim = c(-x_lim, x_lim), ylim = c(-x_lim, x_lim)) + 
  coord_fixed(ratio=1, xlim = c(-x_lim, x_lim), ylim = c(-x_lim, x_lim)) +
  scale_shape_manual(name  ="", values=c(16, 17), labels=c("No or other difference", "AKP down and CHP up"), guide = guide_legend(reverse = TRUE)) +   
  scale_color_manual(name  ="", values=c("black", "red"), labels=c("No or other difference", "AKP down and CHP up"), guide = guide_legend(reverse = TRUE))+
  scale_x_continuous(name = "% AKP Change",
        labels=seq(-25, 25, 5),
        breaks=seq(-x_lim, x_lim, .05)) +
  scale_y_continuous(name = "% CHP Change",
        labels=seq(-25, 25, 5),
        breaks=seq(-x_lim, x_lim, .05)) 

ggsave("figures, appendix/null_1.pdf", width = 10, height = 6, units = "in", dpi = 600)

df_plot %>%
  ggplot(aes(x=AKP_share_diff, y=ABS_share_diff, group=akp_abs_sign)) +
  geom_point(aes(shape=akp_abs_sign, col=akp_abs_sign, shape=akp_abs_sign)) +
  geom_hline(yintercept=0, size=1) +
  geom_vline(xintercept=0, size=1) +
  theme_minimal() + 
  coord_cartesian(xlim = c(-x_lim, x_lim), ylim = c(-x_lim, x_lim)) + 
  coord_fixed(ratio=1, xlim = c(-x_lim, x_lim), ylim = c(-x_lim, x_lim)) +
  scale_shape_manual(name  ="", values=c(16, 17), labels=c("No or other difference", "AKP down and ABS up"), guide = guide_legend(reverse = TRUE)) +   
  scale_color_manual(name  ="", values=c("black", "red"), labels=c("No or other difference", "AKP down and ABS up"), guide = guide_legend(reverse = TRUE))+
  scale_x_continuous(name = "% AKP Change",
        labels=seq(-25, 25, 5),
        breaks=seq(-x_lim, x_lim, .05)) +
  scale_y_continuous(name = "% ABS Change",
        labels=seq(-25, 25, 5),
        breaks=seq(-x_lim, x_lim, .05)) 

ggsave("figures, appendix/null_2.pdf", width = 10, height = 6, units = "in", dpi = 600)

df_plot %>%
  ggplot(aes(x=CHP_share_diff, y=ABS_share_diff, group=chp_abs_sign)) +
  geom_point(aes(shape=chp_abs_sign, col=chp_abs_sign, shape=chp_abs_sign)) +
  geom_hline(yintercept=0, size=1) +
  geom_vline(xintercept=0, size=1) +
  theme_minimal() + 
  coord_cartesian(xlim = c(-x_lim, x_lim), ylim = c(-x_lim, x_lim)) + 
  coord_fixed(ratio=1, xlim = c(-x_lim, x_lim), ylim = c(-x_lim, x_lim)) +
  scale_shape_manual(name  ="", values=c(16, 17), labels=c("No or other difference", "CHP up and ABS down"), guide = guide_legend(reverse = TRUE)) +   
  scale_color_manual(name  ="", values=c("black", "red"), labels=c("No or other difference", "CHP up and ABS down"), guide = guide_legend(reverse = TRUE))+
  scale_x_continuous(name = "% CHP Change",
        labels=seq(-25, 25, 5),
        breaks=seq(-x_lim, x_lim, .05)) +
  scale_y_continuous(name = "% ABS Change",
        labels=seq(-25, 25, 5),
        breaks=seq(-x_lim, x_lim, .05))

ggsave("figures, appendix/null_3.pdf", width = 10, height = 6, units = "in", dpi = 600)
```

# Table 4.1 
* Election results: Heterogeneity in punishment by neighborhood-level covariates

### Prepare the data
```{R}
rm(list=ls(all=TRUE))
set.seed(77)
df_march <- read_csv("data/istanbul, march.csv")
df_june <- read_csv("data/istanbul, june.csv") 
df <- bind_rows(df_march, df_june)

results <- df %>% select(district_name, neighborhood_id, neighborhood_name, box_number, month, CHP, AKP, validvote_total, pop_vote)
rm(df)

# DEFINE ABSTENTION
results <- results %>%  mutate(ABS = pop_vote - (AKP + CHP))

results %>% 
  group_by(district_name, neighborhood_id, neighborhood_name, month) %>% 
  summarize(pop_vote=sum(pop_vote), validvote_total=sum(validvote_total), AKP=sum(AKP), CHP=sum(CHP), ABS=sum(ABS)) -> results_neighbourhood
rm(results)

results_neighbourhood <- results_neighbourhood %>% arrange(neighborhood_id)
results_neighbourhood_march <- results_neighbourhood %>% filter(month=="march")
results_neighbourhood_june <- results_neighbourhood %>% filter(month=="june")

colnames(results_neighbourhood_march) <- paste(colnames(results_neighbourhood_march), "_march", sep = "")
colnames(results_neighbourhood_june) <- paste(colnames(results_neighbourhood_june), "_june", sep = "")

results_neighbourhood_combined <- cbind(results_neighbourhood_march, results_neighbourhood_june)
results <- results_neighbourhood_combined %>% select(district_name_march, neighborhood_id_march, neighborhood_name_march, pop_vote_march, validvote_total_march, AKP_march, CHP_march, ABS_march, pop_vote_june, validvote_total_june, AKP_june, CHP_june, ABS_june)

rm(results_neighbourhood_combined)

results <- results %>% rename(district_name = district_name_march, neighborhood_id=neighborhood_id_march, neighborhood_name = neighborhood_name_march, registered_march=pop_vote_march, valid_march=validvote_total_march, registered_june=pop_vote_june, valid_june=validvote_total_june)

# SAVE FOR LATER
results_neighbourhood_combined <- results
rm(results, results_neighbourhood, results_neighbourhood_march, results_neighbourhood_june)

# MERGE WITH COVARIATE DATA BY DISTRICT AND NEIGHBORHOOD
ids <- read_csv("data/YSK_ids.csv")
results_merged <- left_join(results_neighbourhood_combined, ids, by = c("neighborhood_id"))
results_merged <- results_merged %>% relocate(district_name.y, neighborhood_name.y)
rm(ids, results_neighbourhood_combined)

# MERGE WITH MAHALLEM DATA
covariates <- read_csv("data/mahallem data, to merge.csv")
colnames(covariates)
covariates <- rename(covariates, district_name = "YSK District Name", neighborhood_name="YSK Neighborhood Name", muhtarlik_ID="YSK Neighborhood Code", college_educated="college")
colnames(covariates)

# TRANSFORM TO NUMERIC
covariates$age <- as.numeric(covariates$age)
covariates$home_price <- as.numeric(covariates$home_price)
covariates$density <- as.numeric(covariates$density)
covariates$college_educated <- as.numeric(covariates$college_educated)

# MERGE
covariates_merged <- left_join(results_merged, covariates, by = c("muhtarlik_ID"))
colnames(covariates_merged)
rm(results_merged, covariates)

# PREPARE FOR ANALYSIS
covariates_merged <- covariates_merged %>% select(district_name.x, neighborhood_id, neighborhood_name.x, registered_march, valid_march, AKP_march, CHP_march, ABS_march, registered_june, valid_june, AKP_june, CHP_june, ABS_june, age, density, home_price, college_educated)

results_covariates <- covariates_merged
results_covariates <- results_covariates %>% rename(district_name=district_name.x, neighborhood_name=neighborhood_name.x)

load("data/neighborhoods_transitions.RData")
dim(results_neighbourhood_save)
nrow(results_neighbourhood_save)


results_het <- results_covariates %>% 
  filter(neighborhood_id %in% results_neighbourhood_save$neighborhood_id)
```

### REFORMAT THE DATA TO MIRROR THE EXPERIMENTAL HET ANALYSIS
* This means that March is "control", June is "treatment"
```{R}
# ONLY KEEP THE VARS YOU WILL NEED 
results_het <- results_het %>% 
  select(-registered_march, -registered_june, -valid_march, -valid_june)

# STACK
results_het_march <- results_het %>% 
  select(-ends_with("june")) %>%
  rename(AKP = AKP_march, CHP = CHP_march, ABS = ABS_march) %>%
  mutate(treat=0)

results_het_june <- results_het %>% 
  select(-ends_with("march")) %>%
  rename(AKP = AKP_june, CHP = CHP_june, ABS = ABS_june) %>%
  mutate(treat=1)

results_ols <- bind_rows(results_het_march, results_het_june)
```

### ADD THE SHARE OF AKP/CHP/ABSTAIN OVER REGISTERED VOTERS
```{R}
AKP_share <- results_ols$AKP/(results_ols$AKP + results_ols$CHP + results_ols$ABS)
CHP_share <- results_ols$CHP/(results_ols$AKP + results_ols$CHP + results_ols$ABS)
ABS_share <- results_ols$ABS/(results_ols$AKP + results_ols$CHP + results_ols$ABS)
AKP_two <- results_ols$AKP/(results_ols$AKP + results_ols$CHP)

df_ols <- results_ols %>%
  ungroup() %>%
  select(district_name, neighborhood_id, age, density, home_price, college_educated, treat) %>%
  mutate(AKP_share, CHP_share, ABS_share, AKP_two)
```


## REGRESSION ANALYSIS
```{r}
#AKP VOTE SHARE
df_het_AKP <- df_ols %>%
  mutate(AKP_share_log=log(AKP_share)) %>%
  select(AKP_share_log, treat, age, density, home_price, college_educated, district_name)
est_ols_het_AKP <- lm_robust(AKP_share_log ~ treat*(.-district_name), cluster=district_name, data = df_het_AKP)

#CHP VOTE SHARE
df_het_CHP <- df_ols %>%
  mutate(CHP_share_log=log(CHP_share)) %>%
  select(CHP_share_log, treat, age, density, home_price, college_educated, district_name)
est_ols_het_CHP <- lm_robust(CHP_share_log ~ treat*(.-district_name), cluster=district_name, data = df_het_CHP)

#ABSTENTION RATES
df_het_ABS <- df_ols %>%
  mutate(ABS_share_log=log(ABS_share)) %>%
  select(ABS_share_log, treat, age, density, home_price, college_educated, district_name)
est_ols_het_ABS <- lm_robust(ABS_share_log ~ treat*(.-district_name), cluster=district_name, data = df_het_ABS)

# REGRESSION ANALYSIS: AKP TWO-PARTY VOTE SHARE
df_het_AKP_two <- df_ols %>%
  mutate(AKP_two_log=log(AKP_two)) %>%
  select(AKP_two_log, treat, age, density, home_price, college_educated, district_name)
est_ols_het_AKP_two <- lm_robust(AKP_two_log ~ treat*(.-district_name), cluster=district_name, data = df_het_AKP_two)
```

## CONSTRUCT THE TABLE
```{R}
library(modelsummary)

# AKP TWO-PARTY VOTE SHARE
est_ols_het_two_out <- tidy(est_ols_het_AKP_two)

est_ols_het_two_out_alphas <- est_ols_het_two_out %>%
  slice(c(1:2))

est_ols_het_two_out_betas <- est_ols_het_two_out %>%
  filter(!grepl("treat:",term)) %>% 
  slice(-c(1:2))

est_ols_het_two_out_gammas <- est_ols_het_two_out %>%
  filter(grepl("treat:",term))
est_ols_het_two_out_gammas$term <- est_ols_het_two_out_betas$term


table_coefs_alphas <- est_ols_het_two_out_alphas
table_coefs_alphas$term <- c("Intercept", "$D^-$")
table_coefs_betas <- est_ols_het_two_out_betas
table_coefs_betas$term <- c("Average age", "% college educated", "Population density", "Housing prices")   
table_coefs_gammas <- est_ols_het_two_out_gammas
table_coefs_gammas$term <- c("Average age", "% college educated", "Population density", "Housing prices")   
table_stats <- glance(est_ols_het_AKP_two)[1,c(6,1)]

table_out_alphas_two <- list(
  tidy = table_coefs_alphas,
  glance = table_stats)
class(table_out_alphas_two) <- "modelsummary_list"

table_out_betas_two <- list(
  tidy = table_coefs_betas,
  glance = table_stats)
class(table_out_betas_two) <- "modelsummary_list"

table_out_gammas_two <- list(
  tidy = table_coefs_gammas,
  glance = table_stats)
class(table_out_gammas_two) <- "modelsummary_list"

# AKP VOTE
est_ols_het_1_out <- tidy(est_ols_het_AKP)

est_ols_het_1_out_alphas <- est_ols_het_1_out %>%
  slice(c(1:2))

est_ols_het_1_out_betas <- est_ols_het_1_out %>%
  filter(!grepl("treat:",term)) %>% 
  slice(-c(1:2))

est_ols_het_1_out_gammas <- est_ols_het_1_out %>%
  filter(grepl("treat:",term))
est_ols_het_1_out_gammas$term <- est_ols_het_1_out_betas$term


table_coefs_alphas <- est_ols_het_1_out_alphas
table_coefs_alphas$term <- c("Intercept", "$D^-$")
table_coefs_betas <- est_ols_het_1_out_betas
table_coefs_betas$term <- c("Average age", "% college educated", "Population density", "Housing prices")   
table_coefs_gammas <- est_ols_het_1_out_gammas
table_coefs_gammas$term <- c("Average age", "% college educated", "Population density", "Housing prices")   

table_stats <- glance(est_ols_het_AKP)[1,c(6,1)]

table_out_alphas_1 <- list(
  tidy = table_coefs_alphas,
  glance = table_stats)
class(table_out_alphas_1) <- "modelsummary_list"

table_out_betas_1 <- list(
  tidy = table_coefs_betas,
  glance = table_stats)
class(table_out_betas_1) <- "modelsummary_list"

table_out_gammas_1 <- list(
  tidy = table_coefs_gammas,
  glance = table_stats)
class(table_out_gammas_1) <- "modelsummary_list"

# CHP VOTE
est_ols_het_2_out <- tidy(est_ols_het_CHP)

est_ols_het_2_out_alphas <- est_ols_het_2_out %>%
  slice(c(1:2))

est_ols_het_2_out_betas <- est_ols_het_2_out %>%
  filter(!grepl("treat:",term)) %>% 
  slice(-c(1:2))

est_ols_het_2_out_gammas <- est_ols_het_2_out %>%
  filter(grepl("treat:",term))
est_ols_het_2_out_gammas$term <- est_ols_het_2_out_betas$term


table_coefs_alphas <- est_ols_het_2_out_alphas
table_coefs_alphas$term <- c("Intercept", "$D^-$")
table_coefs_betas <- est_ols_het_2_out_betas
table_coefs_betas$term <- c("Average age", "% college educated", "Population density", "Housing prices")   
table_coefs_gammas <- est_ols_het_2_out_gammas
table_coefs_gammas$term <- c("Average age", "% college educated", "Population density", "Housing prices")   

table_stats <- glance(est_ols_het_CHP)[1,c(6,1)]

table_out_alphas_2 <- list(
  tidy = table_coefs_alphas,
  glance = table_stats)
class(table_out_alphas_2) <- "modelsummary_list"

table_out_betas_2 <- list(
  tidy = table_coefs_betas,
  glance = table_stats)
class(table_out_betas_2) <- "modelsummary_list"

table_out_gammas_2 <- list(
  tidy = table_coefs_gammas,
  glance = table_stats)
class(table_out_gammas_2) <- "modelsummary_list"

# ABSTENTIONS
est_ols_het_0_out <- tidy(est_ols_het_ABS)

est_ols_het_0_out_alphas <- est_ols_het_0_out %>%
  slice(c(1:2))

est_ols_het_0_out_betas <- est_ols_het_0_out %>%
  filter(!grepl("treat:",term)) %>% 
  slice(-c(1:2))

est_ols_het_0_out_gammas <- est_ols_het_0_out %>%
  filter(grepl("treat:",term))
est_ols_het_0_out_gammas$term <- est_ols_het_0_out_betas$term


table_coefs_alphas <- est_ols_het_0_out_alphas
table_coefs_alphas$term <- c("Intercept", "$D^-$")
table_coefs_betas <- est_ols_het_0_out_betas
table_coefs_betas$term <- c("Average age", "% college educated", "Population density", "Housing prices")   
table_coefs_gammas <- est_ols_het_0_out_gammas
table_coefs_gammas$term <- c("Average age", "% college educated", "Population density", "Housing prices")   

table_stats <- glance(est_ols_het_ABS)[1,c(6,1)]

table_out_alphas_0 <- list(
  tidy = table_coefs_alphas,
  glance = table_stats)
class(table_out_alphas_0) <- "modelsummary_list"

table_out_betas_0 <- list(
  tidy = table_coefs_betas,
  glance = table_stats)
class(table_out_betas_0) <- "modelsummary_list"

table_out_gammas_0 <- list(
  tidy = table_coefs_gammas,
  glance = table_stats)
class(table_out_gammas_0) <- "modelsummary_list"

# COMBINE
modelsummary(list(table_out_alphas_1, table_out_betas_1, table_out_gammas_1, table_out_alphas_2, table_out_betas_2, table_out_gammas_2, table_out_alphas_0, table_out_betas_0, table_out_gammas_0, table_out_alphas_two, table_out_betas_two, table_out_gammas_two), 
             stars =  c('*' = .1, '**' = .05,  '***' = .01), 
             escape = FALSE)
```